
library(hydroGOF)
library(clusterSim)
library(sqldf)

data=read.csv("C:/in/myd_merge.csv")
data=subset(data,data$station!=3)

## data$aod=data.Normalization(data$aod,type="n5",normalization="column")
## data$temp=data.Normalization(data$temp,type="n5",normalization="column")
## data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
## data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
## data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")
#data1=data
data1=subset(data,data$pm25!=401.13700)

#model1=lm(pm25~I(aod^2) + I(temp^2) + I(avg_temp^2) + I(avg_rh^2) + I(avg_preci_24^2) + aod*temp + aod*avg_temp + aod*avg_rh + aod*avg_preci_24 + temp*avg_temp + temp*avg_rh + temp*avg_preci_24 + avg_temp*avg_rh + avg_temp*avg_preci_24 + avg_rh*avg_preci_24 + aod + temp + avg_temp + avg_rh + avg_preci_24,data1)



#model1=lm(pm25~aod^2 + temp^2 + avg_temp^2 + avg_rh^2 + avg_preci_24^2 + aod*temp + aod*avg_temp + aod*avg_rh + aod*avg_preci_24 + temp*avg_temp + temp*avg_rh + temp*avg_preci_24 + avg_temp*avg_rh + avg_temp*avg_preci_24 + avg_rh*avg_preci_24 + aod + temp + avg_temp + avg_rh + avg_preci_24,data1)
#model1=lm(pm25~aod^2+temp^2+avg_temp^2+avg_rh^2+avg_preci_24^2+aod*temp+temp*avg_temp+avg_temp*avg_rh+avg_rh*avg_preci_24+avg_preci_24*aod,data1)
model1=lm(pm25~log2(aod)+log2(temp)+log2(avg_temp)+log2(avg_rh)+log2(avg_preci_24),data1)
#model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data1)
#model1=lm(data1$pm25~data1$aod+data1$temp+data1$avg_temp+data1$avg_rh+data1$avg_preci_24)

nSample=nrow(data1)
nPredict= length(model1$coefficients)-1

cutoff=4/(nSample-nPredict-1)
cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)
cookValue=cooks.distance(model1)
data1["cookValue"] = cookValue

data2=subset(data1,data1$cookValue<=cook_thresold)

#model2=lm(pm25~I(aod^2) + I(temp^2) + I(avg_temp^2) + I(avg_rh^2) + I(avg_preci_24^2) + aod*temp + aod*avg_temp + aod*avg_rh + aod*avg_preci_24 + temp*avg_temp + temp*avg_rh + temp*avg_preci_24 + avg_temp*avg_rh + avg_temp*avg_preci_24 + avg_rh*avg_preci_24 + aod + temp + avg_temp + avg_rh + avg_preci_24,data2)


#model2=lm(pm25~aod^2 + temp^2 + avg_temp^2 + avg_rh^2 + avg_preci_24^2 + aod*temp + aod*avg_temp + aod*avg_rh + aod*avg_preci_24 + temp*avg_temp + temp*avg_rh + temp*avg_preci_24 + avg_temp*avg_rh + avg_temp*avg_preci_24 + avg_rh*avg_preci_24 + aod + temp + avg_temp + avg_rh + avg_preci_24,data2)
#model2=lm(pm25~aod^2+temp^2+avg_temp^2+avg_rh^2+avg_preci_24^2+aod*temp+temp*avg_temp+avg_temp*avg_rh+avg_rh*avg_preci_24+avg_preci_24*aod,data2)
model2=lm(pm25~log2(aod)+log2(temp)+log2(avg_temp)+log2(avg_rh)+log2(avg_preci_24),data2)
#model2=lm(lnpm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data2)
#model2=lm(data2$pm25~data2$aod+data2$temp+data2$avg_temp+data2$avg_rh+data2$avg_preci_24)

data2["predict"] = model2$fitted.values



pm25model = data2["predict"]
pm25real= data2["pm25"]
model_re=sum(abs(pm25real-pm25model)/pm25real)/nrow(data2)*100
model_rmse=rmse(pm25real,pm25model)
model_r=cor(pm25real,pm25model)

print("All Daily:")
print(paste("Number of sample: ",nrow(data2)))
print(paste("R2:",round(model_r*model_r,3)))
print(paste("RMSE:",round(model_rmse,3)))
print(paste("RE:",round(model_re,3)))


write.csv(data2,"C:/out/myd_log.csv")

write.csv(model2$coefficients,"C:/out/myd_coeff_4.csv")